intl trade 2

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(readr)
library(readxl)
library(stringr)
library(cowplot)

## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp

library(paletteer)
library(ggrepel)

## Warning: package 'ggrepel' was built under R version 4.3.3

Read data

all_nickel <- read.csv("../data/all_nickel_trade_data/all_nickel_2012.csv", row.names = NULL)
for (year in 2013:2024) {
  new_df = read.csv(paste0("../data/all_nickel_trade_data/all_nickel_", year, ".csv"), row.names = NULL)
  all_nickel <- rbind(all_nickel, new_df)
}

all_nickel <- all_nickel %>%
  select(
    Year = refPeriodId,
    Month = refYear,
    CountryCode = reporterCode,
    CountryName = reporterISO,
    Flow = reporterDesc,
    HSCode = isOriginalClassification,
    Value = fobvalue
  )

Filter data

all_nickel_X <- all_nickel %>%
  filter(Flow == "X") %>%
  group_by(Year, Month, CountryCode, CountryName) %>%
  summarise(Exports = sum(Value)) %>%
  ungroup()

## `summarise()` has grouped output by 'Year', 'Month', 'CountryCode'. You can
## override using the `.groups` argument.

all_nickel_yearly_actual <- read.csv("../data/all_nickel_trade_data/all_nickel_yearly.csv", row.names = NULL) %>%
  select(
    Year = refPeriodId,
    CountryCode = reporterCode,
    CountryName = reporterISO,
    Flow = reporterDesc,
    HSCode = isOriginalClassification,
    Value = fobvalue
  ) %>%
  filter(Flow == "X") %>%
  group_by(Year, CountryCode, CountryName) %>%
  summarise(Exports = sum(Value)) %>%
  ungroup()

## `summarise()` has grouped output by 'Year', 'CountryCode'. You can override
## using the `.groups` argument.

Nickel prices

nickel_price <- read.csv("../data/PNICKUSDM_from2011.csv") %>%
  mutate(Date = date(observation_date)) %>%
  filter(Date >= "2012-01-01") %>%
  mutate(Year = as.numeric(substring(observation_date, 1, 4)))

nickel_price_yearly <- nickel_price %>%
  mutate(Year = substring(observation_date, 1, 4)) %>%
  group_by(Year) %>%
  mutate(Year = as.numeric(Year)) %>%
  summarise(MeanP = mean(PNICKUSDM))

CPI_IDN <- read_xls("../data/API_FP.CPI.TOTL_DS2_en_excel_v2_252315.xls",
                    sheet = "Data", range = "A4:BQ270") %>%
  filter(`Country Name` == "Indonesia") %>%
  select(-`Indicator Name`, -`Indicator Code`) %>%
  pivot_longer(`1960`:`2024`, names_to = "Year", values_to = "CPI") %>%
  mutate(Year = as.numeric(Year)) %>%
  filter(Year > 2010)

USD_to_IDR_conversion <- read.csv("../data/bis_dp_search_export_20251127-202846.csv", skip = 2) %>%
  mutate(Year = 2011:2024) %>%
  select(Year, IDR_USD = OBS_VALUE.Value)

real_P <- CPI_IDN %>%
  left_join(nickel_price, by = c("Year")) %>%
  mutate(RealP = PNICKUSDM / CPI * 100) %>%
  drop_na()

real_P_yearly <- CPI_IDN %>%
  left_join(nickel_price_yearly, by = c("Year")) %>%
  mutate(RealP = MeanP / CPI * 100) %>%
  drop_na()

ggplot(real_P) +
  # annotate("rect", xmin=as.Date("2014-01-01"), xmax=as.Date("2016-12-01"),
  #          ymin=-Inf, ymax=Inf, alpha=0.4, fill="darkgrey") +
  # annotate("rect", xmin=as.Date("2020-01-01"), xmax=as.Date("2024-12-31"),
  #          ymin=-Inf, ymax=Inf, alpha=0.4, fill="darkgrey") +
  geom_line(aes(x = Date, y = RealP/1e3), alpha = 0.4) +
  scale_x_date(
    limits = as.Date(c(NA, "2024-12-31")),
    date_breaks = "1 year",
    date_labels = "%Y"
  ) +
  geom_line(
    data = real_P_yearly %>%
      filter(
        Year >= 2012,
        Year <= 2024
        ) %>%
      mutate(Date = as.Date(paste0(Year, "-06-01"))),
    aes(x = Date, y = RealP/1e3),
    colour = "tomato2", size = 1.1
  ) +
  # geom_point(
  #   data = nickel_price_yearly %>%
  #     filter(
  #       Year >= 2012,
  #       Year <= 2025
  #       ) %>%
  #     mutate(Date = as.Date(paste0(Year, "-06-01"))),
  #   aes(x = Date, y = MeanP/1e3),
  #   colour = "tomato2"
  # ) +
  labs(
    title = "Yearly Mean Price of Nickel from 2012 to 2025",
    subtitle = "Price of Nickel Futures on LME, adjusted by Indonesian CPI (100 in 2010)",
    x = "Date", y = "Price (k USD per Metric Ton)",
    caption = "Source: Federal Reserve Bank of St. Louis"
  ) +
  annotate(geom = "point", x = as.Date("2022-03-01"), y = 20803.46/1e3, colour = "tomato2") +
  annotate(geom = "text", x = as.Date("2022-06-01"), y = 20600/1e3, size = 3, hjust = 0,
           label = "Mar 2022: \nPrice surge") +
  theme_minimal() +
  theme(panel.grid.minor = element_blank())

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

max(real_P$RealP)

## [1] 20803.46

Initial plots

# helper function

k_highest_in_period <- function(year, k) {
  res = all_nickel_yearly_actual %>%
    filter(Year == year) %>%
    slice_max(Exports, n = k) %>%
    pull(CountryCode)
  return(res)
}

k_highest_Q_in_period <- function(year, k) {
  res = nickel_Q %>%
    filter(Year == year) %>%
    slice_max(Quantity, n = k) %>%
    pull(CountryCode)
  return(res)
}

# plotting helpers

plot_ban_rects <- function() {
  list(
    annotate("rect", xmin=2013.9, xmax=2017.1, ymin=-Inf, ymax=Inf, alpha=0.4, fill="darkgrey"),
    annotate("rect", xmin=2019.9, xmax=Inf, ymin=-Inf, ymax=Inf, alpha=0.4, fill="darkgrey")
  )
}

plot_event_lines <- function() {
  list(
    geom_vline(xintercept = 2022, linetype = 2),
    annotate(geom = "text", x = 2019.85, y = 7, hjust = 1, size = 2.5,
             label = "2020: Indonesian export\nban starts again"),
    geom_vline(xintercept = 2020, linetype = 2),
    annotate(geom = "text", x = 2022.15, y = 7.85, hjust = 0, size = 2.5,
             label = "2022: Price surge \nin nickel")
  )
}

all_nickel_X_yearly <- all_nickel_X %>%
  group_by(Year, CountryCode) %>%
  summarise(Exports = sum(Exports)) %>%
  ungroup()

## `summarise()` has grouped output by 'Year'. You can override using the
## `.groups` argument.

all_nickel_X_yearly_to_plot <- all_nickel_X_yearly %>%
  filter(CountryCode != "RUS") %>%
  filter(CountryCode %in% k_highest_in_period(2019, 11)) %>%
  left_join(
    data.frame(
      CountryCode = k_highest_in_period(2019, 11),
      rank = seq_along(k_highest_in_period(2019, 1))
      ),
    by = c("CountryCode")
    )

all_nickel_yearly_actual_to_plot <- all_nickel_yearly_actual %>%
  filter(CountryCode != "RUS") %>%
  filter(CountryCode %in% k_highest_in_period(2019, 11)) %>%
  left_join(
    data.frame(
      CountryCode = k_highest_in_period(2019, 11),
      rank = seq_along(k_highest_in_period(2019, 11))
      ),
    by = c("CountryCode")
    )

k_highest_in_period(2019, 11)

##  [1] "USA" "CAN" "RUS" "IDN" "DEU" "AUS" "GBR" "FIN" "NOR" "ZWE" "PHL"

# plotting

ggplot(all_nickel_X_yearly_to_plot) +
  plot_ban_rects() +
  plot_event_lines() +
  geom_line(aes(x = Year, y = Exports/1e9, group = CountryCode, colour = rank)) +
  geom_text(
    data = all_nickel_X_yearly_to_plot %>%
      filter(Year == 2024) %>%
      select(CountryCode, LastX = Exports, rank),
    aes(x = 2024.2, y = LastX/1e9, label = CountryCode),
    hjust = 0, size = 2.5
    ) +
  scale_colour_gradient2(low = "orange", mid = "mediumseagreen", high = "royalblue", midpoint = 6) +
  geom_line(
    data = all_nickel_X_yearly_to_plot %>% filter(CountryCode == "IDN"),
    aes(x = Year, y = Exports/1e9, group = CountryCode),
    colour = "tomato", linewidth = 1
    ) +
  geom_text(
    data = all_nickel_X_yearly_to_plot %>%
      filter(Year == 2024) %>%
      select(CountryCode, LastX = Exports, rank) %>%
      filter(CountryCode == "IDN"),
    aes(x = 2024.2, y = LastX/1e9, label = CountryCode),
    colour = "tomato", hjust = 0, size = 2.5, fontface = "bold"
    ) +
  scale_x_continuous(limits = c(2012, 2024.5), breaks = 2012:2024) +
  labs(
    title = "Nickel Exports of Top Exporting Countries from 2012 to 2024",
    subtitle = "Top 8 Exporters in Dec 2019, excluding Russia",
    x = "Year", y = "Value of Exports (Billion USD)",
    caption = "Source: UN Comtrade Database"
  ) +
  theme_minimal() +
  theme(legend.position = "None")

ggplot(all_nickel_yearly_actual_to_plot) +
  plot_ban_rects() +
  plot_event_lines() +
  geom_line(aes(x = Year, y = Exports/1e9, group = CountryCode, colour = rank)) +
  geom_text(
    data = all_nickel_yearly_actual_to_plot %>%
      filter(Year == 2024) %>%
      select(CountryCode, LastX = Exports, rank),
    aes(x = 2024.2, y = LastX/1e9, label = CountryCode),
    hjust = 0, size = 2.5
    ) +
  scale_colour_gradient2(low = "orange", mid = "mediumseagreen", high = "royalblue", midpoint = 6) +
  geom_line(
    data = all_nickel_yearly_actual_to_plot %>% filter(CountryCode == "IDN"),
    aes(x = Year, y = Exports/1e9, group = CountryCode),
    colour = "tomato2", linewidth = 1
    ) +
  geom_text(
    data = all_nickel_yearly_actual_to_plot %>%
      filter(Year == 2024) %>%
      select(CountryCode, LastX = Exports, rank) %>%
      filter(CountryCode == "IDN"),
    aes(x = 2024.2, y = LastX/1e9, label = CountryCode),
    colour = "tomato2", hjust = 0, size = 2.5, fontface = "bold"
    ) +
  scale_x_continuous(limits = c(2013, 2024.5), breaks = 2013:2024) +
  labs(
    title = "Nickel Exports of Top Exporting Countries from 2012 to 2024",
    subtitle = "Top 10 Exporters in Dec 2019, excluding Russia",
    x = "Year", y = "Value of Exports (Billion USD)",
    caption = "Source: UN Comtrade Database"
  ) +
  theme_minimal() +
  theme(legend.position = "None")

nickel_Q <- all_nickel_yearly_actual %>%
  left_join(nickel_price_yearly, by = "Year") %>%
  mutate(Quantity = Exports / Year)

to_plot_nickel_Q <- nickel_Q %>%
  filter(CountryCode != "RUS") %>%
  filter(CountryCode %in% k_highest_Q_in_period(2019, 9)) %>%
  left_join(
    data.frame(
      CountryCode = k_highest_Q_in_period(2019, 9),
      rank = seq_along(k_highest_Q_in_period(2019, 9))
      ),
    by = c("CountryCode")
    )

ggplot(to_plot_nickel_Q) +
  list(
    geom_vline(xintercept = 2022, linetype = 2),
    annotate(geom = "text", x = 2022.15, y = 3.9, hjust = 0, size = 2.5,
             label = "2022: Price surge \nin nickel")
  ) +
  geom_line(aes(x = Year, y = Quantity/1e6, group = CountryCode, colour = rank)) +
  geom_text(
    data = to_plot_nickel_Q %>%
      filter(Year == 2024) %>%
      select(CountryCode, LastQ = Quantity, rank),
    aes(x = 2024.2, y = LastQ/1e6, label = CountryCode),
    hjust = 0, size = 2.5
    ) +
  scale_colour_gradient2(low = "orange", mid = "mediumseagreen", high = "royalblue", midpoint = 6) +
  geom_line(
    data = to_plot_nickel_Q %>% filter(CountryCode == "IDN"),
    aes(x = Year, y = Quantity/1e6, group = CountryCode),
    colour = "tomato2", linewidth = 1
    ) +
  geom_text(
    data = to_plot_nickel_Q %>%
      filter(Year == 2024) %>%
      select(CountryCode, LastQ = Quantity, rank) %>%
      filter(CountryCode == "IDN"),
    aes(x = 2024.2, y = LastQ/1e6, label = CountryCode),
    colour = "tomato2", hjust = 0, size = 2.5, fontface = "bold"
    ) +
  scale_x_continuous(limits = c(2013, 2024.5), breaks = 2013:2024) +
  labs(
    title = "Annual Nickel Export Quantity From Top Exporting Countries from 2013 to 2024",
    subtitle = "Top 8 Exporters in 2019, excluding Russia",
    x = "Year", y = "Exports (Million Metric Tonnes)",
    caption = "Source: UN Comtrade Database"
  ) +
  theme_minimal() +
  theme(legend.position = "None")

ggplot(nickel_price) +
  annotate("rect", xmin=as.Date("2014-01-01"), xmax=as.Date("2016-12-01"),
           ymin=-Inf, ymax=Inf, alpha=0.4, fill="darkgrey") +
  annotate("rect", xmin=as.Date("2020-01-01"), xmax=as.Date("2025-01-01"),
           ymin=-Inf, ymax=Inf, alpha=0.4, fill="darkgrey") +
  geom_line(aes(x = Date, y = PNICKUSDM/1e3), colour = "royalblue") +
  scale_x_date(
    limits = as.Date(c(NA, "2025-01-01")),
    date_breaks = "1 year",
    date_labels = "%Y"
  ) +
  geom_line(
    data = nickel_price_yearly %>%
      filter(
        Year >= 2012,
        Year <= 2025
        ) %>%
      mutate(Date = as.Date(paste0(Year, "-05-15"))),
    aes(x = Date, y = MeanP/1e3),
    alpha = 0.5
  ) +
  labs(
    title = "Monthly Price of Nickel from 2012 to 2025 (k USD per Metric Ton)",
    subtitle = "Price of Nickel Futures on LME",
    x = "Date", y = "Price"
  ) +
  annotate(geom = "text", x = as.Date("2022-06-01"), y = 33300/1e3, size = 2, hjust = 0,
           label = "Mar 2022: \nPrice surge") +
  theme_minimal()

## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

ggplot(nickel_price) +
  # annotate("rect", xmin=as.Date("2014-01-01"), xmax=as.Date("2016-12-01"),
  #          ymin=-Inf, ymax=Inf, alpha=0.4, fill="darkgrey") +
  # annotate("rect", xmin=as.Date("2020-01-01"), xmax=as.Date("2024-12-31"),
  #          ymin=-Inf, ymax=Inf, alpha=0.4, fill="darkgrey") +
  geom_line(aes(x = Date, y = PNICKUSDM/1e3), alpha = 0.4) +
  scale_x_date(
    limits = as.Date(c(NA, "2024-12-31")),
    date_breaks = "1 year",
    date_labels = "%Y"
  ) +
  geom_line(
    data = nickel_price_yearly %>%
      filter(
        Year >= 2012,
        Year <= 2024
        ) %>%
      mutate(Date = as.Date(paste0(Year, "-06-01"))),
    aes(x = Date, y = MeanP/1e3),
    colour = "tomato2", size = 1.1
  ) +
  # geom_point(
  #   data = nickel_price_yearly %>%
  #     filter(
  #       Year >= 2012,
  #       Year <= 2025
  #       ) %>%
  #     mutate(Date = as.Date(paste0(Year, "-06-01"))),
  #   aes(x = Date, y = MeanP/1e3),
  #   colour = "tomato2"
  # ) +
  labs(
    title = "Yearly Mean Price of Nickel from 2012 to 2025",
    subtitle = "Price of Nickel Futures on LME",
    x = "Date", y = "Price (k USD per Metric Ton)",
    caption = "Source: "
  ) +
  annotate(geom = "point", x = as.Date("2022-03-01"), y = 33924.18/1e3, colour = "tomato2") +
  annotate(geom = "text", x = as.Date("2022-06-01"), y = 33300/1e3, size = 2, hjust = 0,
           label = "Mar 2022: \nPrice surge") +
  theme_minimal() +
  theme(panel.grid.minor = element_blank())

## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).

Calculating Quantities

nickel_Q <- all_nickel_yearly_actual %>%
  left_join(nickel_price_yearly, by = c("Year")) %>%
  mutate(Quantity = Exports / MeanP)

nickel_Q_to_plot <- nickel_Q %>%
  filter(CountryCode %in% k_highest_in_period(2019, 10)) %>%
  left_join(
    data.frame(
      CountryCode = k_highest_Q_in_period(2019, 10),
      rank = seq_along(k_highest_Q_in_period(2019, 10))
      ),
    by = c("CountryCode")
    )

ggplot(nickel_Q_to_plot, aes(x = Year, y = Quantity)) +
  geom_line(aes(group = CountryCode, colour = rank)) +
  geom_text(
    data = nickel_Q_to_plot %>%
      filter(Year == 2024) %>%
      select(CountryCode, LastQ = Quantity, rank),
    aes(x = 2024.2, y = LastQ, label = CountryCode),
    hjust = 0, size = 2.5
    ) +
  scale_colour_gradient2(low = "orange", mid = "mediumseagreen", high = "royalblue", midpoint = 6) +
  scale_x_continuous(limits = c(2013, 2024.5), breaks = 2013:2024) +
  theme_minimal() +
  theme(legend.position = "None")

ggplot(nickel_Q %>%
         filter(CountryCode == "IDN"), aes(x = MeanP, y = Quantity)) +
  geom_text(aes(label = Year))

nickel_Q_monthly <- all_nickel_X %>%
  left_join(nickel_price %>% mutate(
    Year = as.numeric(substring(observation_date, 1, 4)),
    Month = as.numeric(substring(observation_date, 6, 7))
  ), by = c("Year", "Month")) %>%
  mutate(Quantity = Exports / PNICKUSDM)

nickel_Q_monthly_to_plot <- nickel_Q_monthly %>%
  filter(CountryCode %in% k_highest_in_period(2019, 10)) %>%
  left_join(
    data.frame(
      CountryCode = k_highest_in_period(2019, 10),
      rank = seq_along(k_highest_in_period(2019, 10))
      ),
    by = c("CountryCode")
    )

ggplot(nickel_Q_monthly_to_plot, aes(x = Date, y = Quantity)) +
  geom_line(aes(group = CountryCode, colour = rank)) +
  scale_colour_gradient2(low = "orange", mid = "mediumseagreen", high = "royalblue", midpoint = 6) +
  scale_x_date(
    limits = as.Date(c(NA, "2025-01-01")),
    date_breaks = "1 year",
    date_labels = "%Y"
  ) +
  theme_minimal() +
  theme(legend.position = "None")

Run regression

library(fixest)

## Warning: package 'fixest' was built under R version 4.3.3

reg_df <- all_nickel_yearly_actual %>%
  left_join(nickel_price_yearly, by = "Year") %>%
  mutate(
    Quantity = Exports / MeanP,
    ExportBan = ifelse((((Year >= 2014) & (Year <= 2016)) | (Year >= 2020) & (CountryCode == "IDN")), 1, 0),
    COVID_dummy = ifelse(Year == 2020, 1, 0)
  ) %>%
  mutate(
    Int = ExportBan * MeanP
  )

reg_df <- all_nickel_yearly_actual %>%
  left_join(nickel_price_yearly, by = "Year") %>%
  mutate(
    Quantity = Exports / MeanP,
    ExportBan = ifelse(((Year >= 2014) & (CountryCode == "IDN")), 1, 0),
    COVID_dummy = ifelse(Year == 2020, 1, 0)
  ) %>%
  mutate(
    Int = ExportBan * MeanP
  )

extra_df <- reg_df %>%
  group_by(CountryCode) %>%
  mutate(mean_Q = mean(Quantity)) %>%
  ungroup() %>%
  mutate(translated_Q = Quantity - mean_Q) %>%
  group_by(Year) %>%
  mutate(adj = (ifelse(CountryCode == "IDN", 1, 115)) ** 0.5) %>%
  ungroup() %>%
  group_by(Year, I_IDN = CountryCode == "IDN") %>%
  summarise(mean_Q = sum(translated_Q) / mean(adj)) %>%
  pivot_wider(names_from = I_IDN, values_from = mean_Q) %>%
  mutate(diff = `TRUE` - `FALSE`)

## `summarise()` has grouped output by 'Year'. You can override using the
## `.groups` argument.

ggplot(reg_df, aes(x = Year, y = Quantity)) +
  geom_line(aes(group = CountryCode, alpha = ifelse(CountryCode == "IDN", 1, 0.6))) +
  geom_line(data = reg_df %>% 
              filter(CountryCode != "IDN") %>%
              group_by(Year) %>%
              summarise(Mean_Q = mean(Quantity)), aes(x = Year, y = Mean_Q * 10), colour = "tomato") +
  geom_vline(xintercept = 2014, linetype = 2, 201) +
  theme_minimal()

## Warning: `geom_vline()`: Ignoring `mapping` because `xintercept` was provided.

ggplot(reg_df  %>%
              filter(CountryCode %in% k_highest_in_period(2019, 10)), aes(x = Year, y = Quantity)) +
  geom_line(aes(group = CountryCode, alpha = ifelse(CountryCode == "IDN", 1, 0.6))) +
  geom_line(data = reg_df %>% 
              filter(CountryCode != "IDN") %>%
              filter(CountryCode %in% k_highest_in_period(2019, 10)) %>%
              group_by(Year) %>%
              summarise(Mean_Q = mean(Quantity)), aes(x = Year, y = Mean_Q), colour = "tomato") +
  geom_vline(xintercept = 2014, linetype = 2, 201) +
  theme_minimal()

## Warning: `geom_vline()`: Ignoring `mapping` because `xintercept` was provided.

ggplot(reg_df %>% filter(CountryCode == "IDN"), aes(x = Year, y = log(Quantity))) +
  plot_ban_rects() +
  geom_line() +
  geom_line(data = reg_df %>% 
              filter(CountryCode != "IDN") %>%
              filter(CountryCode %in% k_highest_in_period(2019, 8)) %>%
              group_by(Year) %>%
              summarise(Mean_Q = mean(Quantity)), aes(x = Year, y = log(Mean_Q)), colour = "tomato") +
  geom_vline(xintercept = 2014, linetype = 2, 201) +
  scale_x_continuous(breaks = 2013:2024) +
  theme_minimal()

## Warning: `geom_vline()`: Ignoring `mapping` because `xintercept` was provided.

ggplot(extra_df, aes(x = Year, y = diff)) +
  geom_line() +
  geom_line(data = nickel_price_yearly, aes(x = Year, y = (MeanP - mean(MeanP)) * 20), colour = "royalblue", alpha = 0.7) +
  geom_vline(xintercept = 2014, linetype = 2, 201) +
  theme_minimal()

## Warning: `geom_vline()`: Ignoring `mapping` because `xintercept` was provided.

base_feols <- feols(
  Quantity ~ MeanP + ExportBan + Int + COVID_dummy | CountryCode,
  data = reg_df %>% filter(Year > 2013))

## The variable 'ExportBan' has been removed because of collinearity (see $collin.var).

base_feols

## OLS estimation, Dep. Var.: Quantity
## Observations: 1,312
## Fixed-effects: CountryCode: 169
## Standard-errors: Clustered (CountryCode) 
##                 Estimate  Std. Error   t value  Pr(>|t|)    
## MeanP          -0.277131    0.108618  -2.55142  0.011621 *  
## Int            14.220362    0.102693 138.47419 < 2.2e-16 ***
## COVID_dummy -2467.542032 1198.239091  -2.05931  0.041008 *  
## ... 1 variable was removed because of collinearity (ExportBan)
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 15,020.6     Adj. R2: 0.866329
##                  Within R2: 0.143057

base_feols_log_1 <- feols(
  log(Quantity) ~ log(MeanP) + ExportBan + ExportBan * log(MeanP) + COVID_dummy | CountryCode,
  data = reg_df %>% filter(CountryCode %in% k_highest_in_period(2019, 1)))

## The variables 'ExportBan' and 'log(MeanP):ExportBan' have been removed because of collinearity (see $collin.var).

log_1_ts = base_feols_log_1$coefficients / base_feols_log_1$se
df_temp <- data.frame(k = 1, t_p = log_1_ts["log(MeanP)"], t_ban = log_1_ts["ExportBan"], t_int = log_1_ts["log(MeanP):ExportBan"])
for (i in 2:50) {
  feols = feols(
    log(Quantity) ~ log(MeanP) + ExportBan + ExportBan * log(MeanP) + COVID_dummy | CountryCode,
    data = reg_df %>% filter(CountryCode %in% k_highest_in_period(2019, i)))
  Ts = feols$coefficients / feols$se
  new_row = data.frame(
    k = i,
    t_p = Ts["log(MeanP)"],
    t_ban = Ts["ExportBan"],
    t_int = Ts["log(MeanP):ExportBan"]
  )
  df_temp <- rbind(df_temp, new_row)
}

## The variables 'ExportBan' and 'log(MeanP):ExportBan' have been removed because of collinearity (see $collin.var).

## The variables 'ExportBan' and 'log(MeanP):ExportBan' have been removed because of collinearity (see $collin.var).

ggplot(df_temp, aes(x = k)) +
  geom_hline(yintercept = -1.96, colour = "tomato", linetype = 2) +
  geom_hline(yintercept = 1.96, colour = "tomato", linetype = 2) +
  geom_line(aes(y = t_p), colour = "goldenrod") +
  geom_line(aes(y = t_ban), colour = "seagreen") +
  geom_line(aes(y = t_int), colour = "royalblue")

## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 3 rows containing missing values or values outside the scale range
## (`geom_line()`).

Simple regression (only Indonesia)

reg_df_1 <- all_nickel_yearly_actual %>%
  left_join(nickel_price_yearly, by = "Year") %>%
  mutate(
    Quantity = Exports / MeanP,
    After = ifelse(((Year >= 2014) & (Year <= 2016)) | (Year >= 2020), 1, 0),
    ExportBan = ifelse(CountryCode == "IDN", 1, 0),
    COVID_dummy = ifelse(Year == 2020, 1, 0)
  ) %>%
  mutate(
    B1 = After,
    B2 = After * ExportBan,
    B3 = log(MeanP),
    B4 = log(MeanP) * After,
    B5 = log(MeanP) * ExportBan,
    B6 = log(MeanP) * ExportBan * After,
    Y = log(Quantity)
  )


reg_df_1 <- all_nickel_yearly_actual %>%
  left_join(nickel_price_yearly, by = "Year") %>%
  filter(Year >= 2013) %>%
  mutate(
    Quantity = Exports / MeanP,
    After = ifelse(Year >= 2020, 1, 0),
    Indo = ifelse(CountryCode == "IDN", 1, 0),
    COVID_dummy = ifelse(Year == 2020, 1, 0)
  ) %>%
  mutate(
    B1 = After,
    B2 = Indo,
    B3 = After * Indo,
    B4 = log(MeanP),
    B5 = log(MeanP) * After,
    B6 = log(MeanP) * Indo,
    B7 = log(MeanP) * Indo * After,
    Y = log(Quantity)
  )


base_feols_log <- feols(
  Y ~ B1 + B3 + B4 + B6 + B7 | CountryCode,
  data = reg_df_1)
base_feols_log

## OLS estimation, Dep. Var.: Y
## Observations: 1,428
## Fixed-effects: CountryCode: 172
## Standard-errors: Clustered (CountryCode) 
##     Estimate   Std. Error       t value   Pr(>|t|)    
## B1  0.084486 1.770513e-01  4.771850e-01 6.3384e-01    
## B3 -8.228401 1.770513e-01 -4.647466e+01  < 2.2e-16 ***
## B4 -0.069190 1.857459e-01 -3.725000e-01 7.0998e-01    
## B6  1.014670 1.857459e-01  5.462678e+00 1.6322e-07 ***
## B7  0.851569 2.840000e-14  2.999195e+13  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 1.61691     Adj. R2: 0.875782
##                 Within R2: 8.945e-4

base_feols_log <- feols(
  Y ~ B1 + B3 + B4 + B5 + B6 + B7 | CountryCode,
  data = reg_df_1 %>% filter(CountryCode != "RUS")
  # %>% filter(CountryCode %in% k_highest_in_period(2019, 10))
  , vcov = "hc1")
etable(base_feols_log)

##                    base_feols_log
## Dependent Var.:                 Y
##                                  
## B1               -12.43** (4.326)
## B3                  4.290 (14.07)
## B4              -0.6902* (0.3266)
## B5               1.296** (0.4484)
## B6                1.636. (0.8875)
## B7                -0.4448 (1.440)
## Fixed-Effects:  -----------------
## CountryCode                   Yes
## _______________ _________________
## S.E. type       Heteroskeda.-rob.
## Observations                1,419
## R2                        0.89005
## Within R2                 0.00690
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# robust to COVID control

base_feols_log <- feols(
  Y ~ B1 + B3 + B4 + B5 + B6 + B7,
  data = reg_df_1 %>% filter(CountryCode != "RUS")
  , vcov = "hc1")
etable(base_feols_log)

##                     base_feols_log
## Dependent Var.:                  Y
##                                   
## Constant             6.559 (8.473)
## B1                  -6.678 (12.38)
## B3                  -5.531 (13.67)
## B4                -0.2337 (0.8970)
## B5                  0.7109 (1.281)
## B6              0.7494*** (0.0226)
## B7                  0.5705 (1.373)
## _______________ __________________
## S.E. type       Heteroskedas.-rob.
## Observations                 1,419
## R2                         0.01931
## Adj. R2                    0.01515
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

base_feols_log <- feols(
  Y ~ B3 + B4 | CountryCode + Year,
  data = reg_df_1 %>% filter(!CountryCode %in% c("RUS", "PHL"))
  , vcov = "hc1")

## The variable 'B4' has been removed because of collinearity (see $collin.var).

base_feols_log <- feols(
  Y ~ B3 | CountryCode + Year,
  data = reg_df_1 %>% filter(!CountryCode %in% c("RUS")))
etable(base_feols_log)

##                     base_feols_log
## Dependent Var.:                  Y
##                                   
## B3              0.5419*** (0.1602)
## Fixed-Effects:  ------------------
## CountryCode                    Yes
## Year                           Yes
## _______________ __________________
## S.E.: Clustered    by: CountryCode
## Observations                 1,419
## R2                         0.89085
## Within R2                  0.00023
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

base_feols_log <- feols(
  Y ~ B3 + B7 | CountryCode + Year,
  data = reg_df_1 %>% filter(!CountryCode %in% c("RUS")))
etable(base_feols_log)

##                    base_feols_log
## Dependent Var.:                 Y
##                                  
## B3              -11.34*** (2.901)
## B7              1.207*** (0.2949)
## Fixed-Effects:  -----------------
## CountryCode                   Yes
## Year                          Yes
## _______________ _________________
## S.E.: Clustered   by: CountryCode
## Observations                1,419
## R2                        0.89086
## Within R2                 0.00032
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lm1 <- lm(Quantity ~ MeanP + After + I(MeanP * After) + COVID_dummy, data = reg_df_1 %>%
            filter(CountryCode == "IDN") %>%
            mutate(After = Year >= 2022))
summary(lm1)

## 
## Call:
## lm(formula = Quantity ~ MeanP + After + I(MeanP * After) + COVID_dummy, 
##     data = reg_df_1 %>% filter(CountryCode == "IDN") %>% mutate(After = Year >= 
##         2022))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -36422 -23966  -5748  11498  76238 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       65978.519  71617.087   0.921 0.387559    
## MeanP                 2.166      5.130   0.422 0.685539    
## AfterTRUE        858224.099 158620.130   5.411 0.000997 ***
## I(MeanP * After)    -29.432      8.293  -3.549 0.009357 ** 
## COVID_dummy      -37226.724  44251.115  -0.841 0.428001    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 41720 on 7 degrees of freedom
## Multiple R-squared:  0.9337, Adjusted R-squared:  0.8959 
## F-statistic: 24.66 on 4 and 7 DF,  p-value: 0.0003196

lm1 <- lm(log(Quantity) ~ log(MeanP) + After + I(log(MeanP) * After) + COVID_dummy, data = reg_df_1 %>%
            filter(CountryCode == "IDN") %>%
            mutate(After = Year >= 2022))
summary(lm1)

## 
## Call:
## lm(formula = log(Quantity) ~ log(MeanP) + After + I(log(MeanP) * 
##     After) + COVID_dummy, data = reg_df_1 %>% filter(CountryCode == 
##     "IDN") %>% mutate(After = Year >= 2022))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.37025 -0.22855 -0.00354  0.05554  0.62909 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)
## (Intercept)             7.7996     5.9376   1.314    0.230
## log(MeanP)              0.3787     0.6249   0.606    0.564
## AfterTRUE              21.7010    13.6080   1.595    0.155
## I(log(MeanP) * After)  -2.0668     1.3793  -1.498    0.178
## COVID_dummy            -0.4307     0.3991  -1.079    0.316
## 
## Residual standard error: 0.3758 on 7 degrees of freedom
## Multiple R-squared:  0.8209, Adjusted R-squared:  0.7186 
## F-statistic: 8.023 on 4 and 7 DF,  p-value: 0.009411

base_feols_log <- feols(
  log(Quantity) ~ log(MeanP) + ExportBan + ExportBan * log(MeanP) + COVID_dummy,
  data = reg_df %>% filter(CountryCode %in% k_highest_in_period(2019, 10)),
  vcov = "hc1")
base_feols_log

## OLS estimation, Dep. Var.: log(Quantity)
## Observations: 117
## Standard-errors: Heteroskedasticity-robust 
##                        Estimate Std. Error   t value   Pr(>|t|)    
## (Intercept)           15.764421   2.766217  5.698909 9.8986e-08 ***
## log(MeanP)            -0.436004   0.289096 -1.508163 1.3433e-01    
## ExportBan            -18.079726   4.225375 -4.278845 3.9758e-05 ***
## COVID_dummy           -0.030677   0.178009 -0.172335 8.6348e-01    
## log(MeanP):ExportBan   1.888710   0.450143  4.195802 5.4666e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.648308   Adj. R2: 0.032696

base_feols_log$coefficients

##          (Intercept)           log(MeanP)            ExportBan 
##          15.76442127          -0.43600404         -18.07972613 
##          COVID_dummy log(MeanP):ExportBan 
##          -0.03067732           1.88870974

base_feols_log$se

##          (Intercept)           log(MeanP)            ExportBan 
##            2.7662174            0.2890960            4.2253751 
##          COVID_dummy log(MeanP):ExportBan 
##            0.1780094            0.4501427 
## attr(,"type")
## [1] "Heteroskedasticity-robust"

Reg formula:

\[ ln(Q) = \beta_{0} + \beta_{1}After_{t} + \beta_{2}After_{t} * ExportBan_{i} + \beta_{3}ln(P_{t}) + \beta_{4}ln(P_{t}) * After_{t} + \beta_{5}ln(Price_{t}) * ExportBan_{i} + \beta_{6}(all 3) \]

reg_df_3 <- all_nickel_yearly_actual %>%
  left_join(real_P_yearly, by = "Year") %>%
  filter(Year >= 2013) %>%
  mutate(
    Quantity = Exports / RealP,
    After = ifelse(Year >= 2022, 1, 0),
    Indo = ifelse(CountryCode == "IDN", 1, 0),
    COVID_dummy = ifelse(Year == 2020, 1, 0)
  ) %>%
  mutate(
    B1 = After,
    B2 = Indo,
    B3 = After * Indo,
    B4 = log(RealP),
    B5 = log(RealP) * After,
    B6 = log(RealP) * Indo,
    B7 = log(RealP) * Indo * After,
    Y = log(Quantity)
  )

feols_log_3 <- feols(
  Y ~ B4 + B5 + B6 + B7 | CountryCode,
  data = reg_df_3 %>% filter(!CountryCode %in% c("RUS")))
etable(feols_log_3)

##                        feols_log_3
## Dependent Var.:                  Y
##                                   
## B4              -0.6432** (0.2125)
## B5              0.0666*** (0.0191)
## B6                0.5283* (0.2125)
## B7              0.0978*** (0.0191)
## Fixed-Effects:  ------------------
## CountryCode                    Yes
## _______________ __________________
## S.E.: Clustered    by: CountryCode
## Observations                 1,419
## R2                         0.89042
## Within R2                  0.02079
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

feols_log_event_study_0 <- feols(
  Y ~ i(Year, ref = 2021) * Indo | CountryCode,
  data = reg_df_3 %>% filter(!CountryCode %in% c("RUS")))

## The variable 'Indo' has been removed because of collinearity (see $collin.var).

etable(feols_log_event_study_0)

##                    feols_log_event..0
## Dependent Var.:                     Y
##                                      
## Year = 2013          -0.3945 (0.2970)
## Year = 2014          -0.2402 (0.2902)
## Year = 2015          -0.2361 (0.2820)
## Year = 2016           0.2755 (0.2358)
## Year = 2017           0.1668 (0.2316)
## Year = 2018           0.0020 (0.1968)
## Year = 2019           0.1177 (0.2069)
## Year = 2020          -0.1629 (0.1827)
## Year = 2022           0.2480 (0.1574)
## Year = 2023         0.5366** (0.1880)
## Year = 2024          0.4216. (0.2491)
## Indo x Year = 2013  1.024*** (0.2970)
## Indo x Year = 2014   -0.0161 (0.2902)
## Indo x Year = 2015    0.0449 (0.2820)
## Indo x Year = 2016  -0.5227* (0.2358)
## Indo x Year = 2017   -0.1606 (0.2316)
## Indo x Year = 2018   0.3752. (0.1968)
## Indo x Year = 2019   0.5277* (0.2069)
## Indo x Year = 2020   -0.0237 (0.1827)
## Indo x Year = 2022 0.9876*** (0.1574)
## Indo x Year = 2023  1.056*** (0.1880)
## Indo x Year = 2024  1.600*** (0.2491)
## Fixed-Effects:     ------------------
## CountryCode                       Yes
## __________________ __________________
## S.E.: Clustered       by: CountryCode
## Observations                    1,419
## R2                            0.89116
## Within R2                     0.02740
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

feols_log_event_study <- feols(
  Y ~ i(Year, ref = 2021) * Indo | CountryCode + Year,
  data = reg_df_3 %>% filter(!CountryCode %in% c("RUS")))

## The variables 'Year::2013', 'Year::2014', 'Year::2015', 'Year::2016', 'Year::2017', 'Year::2018' and 6 others have been removed because of collinearity (see $collin.var).

etable(feols_log_event_study)

##                    feols_log_event_..
## Dependent Var.:                     Y
##                                      
## Indo x Year = 2013  1.024*** (0.2970)
## Indo x Year = 2014   -0.0161 (0.2902)
## Indo x Year = 2015    0.0449 (0.2820)
## Indo x Year = 2016  -0.5227* (0.2358)
## Indo x Year = 2017   -0.1606 (0.2316)
## Indo x Year = 2018   0.3752. (0.1968)
## Indo x Year = 2019   0.5277* (0.2069)
## Indo x Year = 2020   -0.0237 (0.1827)
## Indo x Year = 2022 0.9876*** (0.1574)
## Indo x Year = 2023  1.056*** (0.1880)
## Indo x Year = 2024  1.600*** (0.2491)
## Fixed-Effects:     ------------------
## CountryCode                       Yes
## Year                              Yes
## __________________ __________________
## S.E.: Clustered       by: CountryCode
## Observations                    1,419
## R2                            0.89116
## Within R2                     0.00119
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

event_study_plotting_data <- data.frame(
  Coefficient = feols_log_event_study$coefficients,
  se = feols_log_event_study_0$se[1:11]
) %>%
  mutate(Year = c(2013:2020, 2022:2024)) %>%
  mutate(p95_length = se * 1.96)

ref_row <- data.frame(
  Year       = 2021,
  Coefficient = 0,
  se = 0,
  p95_length = 0
)
plot_df <- rbind(event_study_plotting_data, ref_row)

ggplot(plot_df, aes(
  x = Year,
  y = Coefficient,
  ymin = Coefficient - p95_length,
  ymax = Coefficient + p95_length
  )) + 
  geom_hline(yintercept = 0, colour = "grey") +
  geom_vline(xintercept = 2021, linetype = 2) +
  geom_pointrange() +
  geom_line() +
  scale_x_continuous(breaks = 2013:2024) +
  labs(title = "Event Study of Effect of Downstreaming Policies on Export Level") +
  theme_minimal()

feols_log_event_study_2 <- feols(
  Y ~ i(Year, ref = 2021) * Indo * log(MeanP) | CountryCode,
  data = reg_df_1 %>% filter(!CountryCode %in% c("RUS")))

## The variables 'Indo', 'log(MeanP)', 'log(MeanP):Year::2013', 'log(MeanP):Year::2015', 'log(MeanP):Year::2016', 'log(MeanP):Year::2017' and 18 others have been removed because of collinearity (see $collin.var).

etable(feols_log_event_study_2)

##                          feols_log_event..2
## Dependent Var.:                           Y
##                                            
## Year = 2013                -0.1029 (0.2972)
## Year = 2014                -0.0104 (0.2518)
## Year = 2015                -0.0683 (0.2821)
## Year = 2016                0.4087. (0.2359)
## Year = 2017                 0.2627 (0.2316)
## Year = 2018                 0.0664 (0.1968)
## Year = 2019                 0.1522 (0.2070)
## Year = 2020                -0.1474 (0.1827)
## Year = 2022                 0.2068 (0.1574)
## Year = 2023                0.4593* (0.1881)
## Year = 2024                 0.3228 (0.2492)
## Indo x Year = 2013        1.024*** (0.2972)
## Indo x Year = 2014         -0.0161 (0.2903)
## Indo x Year = 2015          0.0449 (0.2821)
## Indo x Year = 2016        -0.5227* (0.2359)
## Indo x Year = 2017         -0.1606 (0.2316)
## Indo x Year = 2018         0.3752. (0.1968)
## Indo x Year = 2019         0.5277* (0.2070)
## Indo x Year = 2020         -0.0237 (0.1827)
## Indo x Year = 2022       0.9876*** (0.1574)
## Indo x Year = 2023        1.056*** (0.1881)
## Indo x Year = 2024        1.600*** (0.2492)
## log(MeanP) x Year = 2014  -3.24e-5 (0.0045)
## Fixed-Effects:           ------------------
## CountryCode                             Yes
## ________________________ __________________
## S.E.: Clustered             by: CountryCode
## Observations                          1,419
## R2                                  0.89096
## Within R2                           0.01509
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

event_study_plotting_data <- data.frame(
  Coefficient = feols_log_event_study$coefficients[12:22],
  se = feols_log_event_study$se[1:11]
) %>%
  mutate(Year = c(2013:2020, 2022:2024)) %>%
  mutate(p95_length = se * 1.96)

ref_row <- data.frame(
  Year       = 2021,
  Coefficient = 0,
  se = 0,
  p95_length = 0
)
plot_df <- rbind(event_study_plotting_data, ref_row)

ggplot(plot_df, aes(
  x = Year,
  y = Coefficient,
  ymin = Coefficient - p95_length,
  ymax = Coefficient + p95_length
  )) + 
  geom_hline(yintercept = 0, colour = "grey") +
  geom_vline(xintercept = 2021, linetype = 2) +
  geom_pointrange() +
  geom_line() +
  scale_x_continuous(breaks = 2013:2024) +
  theme_minimal()

## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_pointrange()`).

## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_line()`).

## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

reg_df_2 <- all_nickel_yearly_actual %>%
  left_join(real_P_yearly, by = "Year") %>%
  filter(Year >= 2013) %>%
  group_by(CountryCode) %>%
  mutate(
    Qt = Exports / RealP,
    Pt = RealP,
  ) %>%
  mutate(
    Qt_1 = lag(Qt),
    Pt_1 = lag(Pt)
  ) %>%
  ungroup() %>%
  mutate(
    Indo = ifelse(CountryCode == "IDN", 1, 0),
    Downstream = ifelse(Year >= 2020, 1, 0),
    After = ifelse(Year >= 2022, 1, 0),
    COVID_dummy = ifelse(Year == 2020, 1, 0),
    Y = log(Qt / Qt_1),
    dP = log(Pt / Pt_1)
  ) %>%
  drop_na()

feols_log_2 <- feols(
  Y ~ Indo * Downstream * dP | CountryCode,
  data = reg_df_2 %>% filter(!CountryCode %in% c("RUS")) %>% filter(CountryCode %in% k_highest_in_period(2019, 20)))

## The variable 'Indo' has been removed because of collinearity (see $collin.var).

etable(feols_log_2)

##                               feols_log_2
## Dependent Var.:                         Y
##                                          
## Downstream              -0.1341. (0.0696)
## dP                        0.6848 (0.7234)
## Indo x Downstream      0.3907*** (0.0696)
## Indo x dP                -0.5617 (0.7234)
## Downstream x dP          -0.9541 (0.7666)
## Indo x Downstream x dP    1.644* (0.7666)
## Fixed-Effects:         ------------------
## CountryCode                           Yes
## ______________________ __________________
## S.E.: Clustered           by: CountryCode
## Observations                          201
## R2                                0.00718
## Within R2                         0.00591
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

wald(feols_log_2, "dP + dP:Indo = 0")

## [1] NA

Number of smelters

https://pubs.usgs.gov/myb/vol3/2023/myb3-2023-indonesia.pdf

num_smelters <- data.frame(
  Year = c(2013, 2014, 2020, 2023),
  completed = c(2, 3, 13, 44),
  constructing = c(0, 0, 0, 44+28),
  planned = c(0, 0, 0, 44+28+24)
)

num_complete <- num_smelters %>% select(Year, completed)
for (i in 2015:2024) {
  if (i %in% num_smelters$Year) {next}
  else if (i < 2020) {
    num_complete <- rbind(
      num_complete,
      data.frame(
        Year = i,
        completed = num_smelters[2, "completed"] + (num_smelters[3, "completed"] - num_smelters[2, "completed"]) * (i - 2014) / 6
      )
    )
  }
  else if (i < 2023) {
    num_complete <- rbind(
      num_complete,
      data.frame(
        Year = i,
        completed = num_smelters[3, "completed"] + (num_smelters[4, "completed"] - num_smelters[3, "completed"]) * (i - 2020) / 3
      )
    )
  }
  else {
    num_complete <- rbind(
      num_complete,
      data.frame(
        Year = i,
        completed = num_smelters[4, "completed"] + (num_smelters[4, "completed"] - num_smelters[3, "completed"]) * (i - 2023) / 3
      )
    )
  }
}

ggplot(num_smelters, aes(x = Year)) +
  geom_col(aes(y = completed), fill = "royalblue") +
  geom_point(aes(y = completed)) +
  geom_line(aes(y = completed)) +
  geom_text(aes(y = completed + 3, label = completed)) +
  labs(title = "Number of Nickel Smelters in Indonesia",
       x = "Year", y = "Number of Smelters",
       caption = "Source: Chung, 2023") +
  scale_x_continuous(breaks = 2013:2024) +
  theme_minimal()

ggplot(num_complete, aes(x = Year)) +
  geom_col(aes(y = completed), fill = "royalblue") +
  geom_point(aes(y = completed)) +
  geom_line(aes(y = completed)) +
  geom_text(aes(y = completed + 3, label = completed)) +
  labs(title = "Number of Nickel Smelters in Indonesia",
       x = "Year", y = "Number of Smelters",
       caption = "Source: ") +
  scale_x_continuous(breaks = 2013:2024) +
  theme_minimal()

reg_df_4 <- all_nickel_yearly_actual %>%
  left_join(real_P_yearly, by = "Year") %>%
  filter(Year >= 2014) %>%
  mutate(
    Quantity = Exports / RealP,
    Downstream = ifelse(Year >= 2020, 1, 0),
    Surge = ifelse(Year %in% c(2022), 1, 0),
    Decline = ifelse(Year %in% c(2023, 2024), 1, 0),
    After = ifelse(Year >= 2022, 1, 0),
    Indo = ifelse(CountryCode == "IDN", 1, 0),
    COVID_dummy = ifelse(Year == 2020, 1, 0)
  ) %>%
  left_join(num_complete, by = "Year")

feols_log_4 <- feols(
  log(Quantity) ~ Surge * log(RealP) + Decline * log(RealP) + Indo * Downstream | CountryCode,
  data = reg_df_4
)

## The variables 'Indo' and 'Surge:log(RealP)' have been removed because of collinearity (see $collin.var).

etable(feols_log_4)

##                             feols_log_4
## Dependent Var.:           log(Quantity)
##                                        
## Surge                 0.5863** (0.1785)
## log(RealP)            -0.5601* (0.2550)
## Decline                  -8.067 (7.006)
## Downstream             -0.0347 (0.1664)
## log(RealP) x Decline    0.9324 (0.7459)
## Indo x Downstream    0.6899*** (0.1503)
## Fixed-Effects:       ------------------
## CountryCode                         Yes
## ____________________ __________________
## S.E.: Clustered         by: CountryCode
## Observations                      1,312
## R2                              0.89721
## Within R2                       0.01864
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

feols_log_4 <- feols(
  log(Quantity) ~ Indo * After * log(RealP) + Indo * Downstream | CountryCode,
  data = reg_df_4
)

## The variable 'Indo' has been removed because of collinearity (see $collin.var).

etable(feols_log_4)

##                                 feols_log_4
## Dependent Var.:               log(Quantity)
##                                            
## After                        -1.408 (5.411)
## log(RealP)                -0.5644* (0.2572)
## Downstream                 -0.0277 (0.1668)
## Indo x After               18.16*** (5.411)
## Indo x log(RealP)          0.5496* (0.2572)
## After x log(RealP)          0.2130 (0.5705)
## Indo x Downstream          -0.1192 (0.1668)
## Indo x After x log(RealP) -1.807** (0.5705)
## Fixed-Effects:            -----------------
## CountryCode                             Yes
## _________________________ _________________
## S.E.: Clustered             by: CountryCode
## Observations                          1,312
## R2                                  0.89718
## Within R2                           0.01834
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

feols_log_4 <- feols(
  log(Quantity) ~ After * log(RealP) + completed,
  data = reg_df_4 %>% filter(CountryCode == "IDN")
)
etable(feols_log_4)

##                         feols_log_4
## Dependent Var.:       log(Quantity)
##                                    
## Constant             13.27* (5.377)
## After                 8.402 (14.64)
## log(RealP)         -0.1980 (0.5929)
## completed           0.0152 (0.0206)
## After x log(RealP)  -0.7696 (1.511)
## __________________ ________________
## S.E. type                       IID
## Observations                     11
## R2                          0.89631
## Adj. R2                     0.82718
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

feols_log_4 <- feols(
  log(Quantity) ~ After * log(RealP) + completed * Indo + COVID_dummy | CountryCode,
  data = reg_df_4
)

## The variable 'Indo' has been removed because of collinearity (see $collin.var).

etable(feols_log_4)

##                           feols_log_4
## Dependent Var.:         log(Quantity)
##                                      
## After                  -7.554 (9.188)
## log(RealP)          -0.6719* (0.2638)
## completed             0.0117 (0.0140)
## COVID_dummy         -0.2516. (0.1410)
## After x log(RealP)    0.8202 (0.9386)
## completed x Indo   0.0321*** (0.0048)
## Fixed-Effects:     ------------------
## CountryCode                       Yes
## __________________ __________________
## S.E.: Clustered       by: CountryCode
## Observations                    1,312
## R2                            0.89744
## Within R2                     0.02087
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

feols_log_4 <- feols(
  log(Quantity) ~ log(RealP) + Indo * Downstream | CountryCode,
  data = reg_df_4
)

## The variable 'Indo' has been removed because of collinearity (see $collin.var).

etable(feols_log_4)

##                          feols_log_4
## Dependent Var.:        log(Quantity)
##                                     
## log(RealP)          -0.2030 (0.1972)
## Downstream           0.2311 (0.1657)
## Indo x Downstream 0.6981*** (0.1493)
## Fixed-Effects:    ------------------
## CountryCode                      Yes
## _________________ __________________
## S.E.: Clustered      by: CountryCode
## Observations                   1,312
## R2                           0.89570
## Within R2                    0.00425
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# event study

feols_log_4_event_study <- feols(
  log(Quantity) ~ log(RealP) * i(Year, ref = 2021) + Indo * i(Year, ref = 2021) | CountryCode,
  data = reg_df_4
)

## The variables 'Year::2024', 'Indo', 'log(RealP):Year::2016', 'log(RealP):Year::2018', 'log(RealP):Year::2020', 'log(RealP):Year::2022' and 2 others have been removed because of collinearity (see $collin.var).

etable(feols_log_4_event_study)

##                          feols_log_4_even..
## Dependent Var.:               log(Quantity)
##                                            
## log(RealP)                  -2.230. (1.279)
## Year = 2014                 0.2318 (0.2989)
## Year = 2015               -0.5113* (0.2210)
## Year = 2016                -0.8807 (0.6358)
## Year = 2017                -0.4850 (0.3922)
## Year = 2018               -0.6164. (0.3578)
## Year = 2019                -0.2812 (0.2000)
## Year = 2020               -0.7795* (0.3288)
## Year = 2022                0.9138. (0.4637)
## Year = 2023               0.7071** (0.2608)
## log(RealP) x Year = 2014   -0.0158 (0.0121)
## log(RealP) x Year = 2015  -0.0351* (0.0150)
## log(RealP) x Year = 2017   -0.0452 (0.0276)
## log(RealP) x Year = 2019   -0.0161 (0.0116)
## Indo x Year = 2014         -0.0242 (0.2891)
## Indo x Year = 2015          0.0288 (0.2804)
## Indo x Year = 2016        -0.5293* (0.2318)
## Indo x Year = 2017         -0.1714 (0.2267)
## Indo x Year = 2018         0.3739. (0.1904)
## Indo x Year = 2019         0.5193* (0.2029)
## Indo x Year = 2020         -0.0236 (0.1814)
## Indo x Year = 2022       0.9813*** (0.1562)
## Indo x Year = 2023        1.055*** (0.1848)
## Indo x Year = 2024        1.592*** (0.2463)
## Fixed-Effects:           ------------------
## CountryCode                             Yes
## ________________________ __________________
## S.E.: Clustered             by: CountryCode
## Observations                          1,312
## R2                                  0.89782
## Within R2                           0.02444
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

event_study_plotting_data <- data.frame(
  Coefficient = feols_log_4_event_study$coefficients[12:22],
  se = feols_log_4_event_study$se[1:11]
) %>%
  mutate(Year = c(2013:2020, 2022:2024)) %>%
  mutate(p95_length = se * 1.96)

ref_row <- data.frame(
  Year       = 2021,
  Coefficient = 0,
  se = 0,
  p95_length = 0
)
plot_df <- rbind(event_study_plotting_data, ref_row)

ggplot(plot_df, aes(
  x = Year,
  y = Coefficient,
  ymin = Coefficient - p95_length,
  ymax = Coefficient + p95_length
  )) + 
  geom_hline(yintercept = 0, colour = "grey") +
  geom_vline(xintercept = 2021, linetype = 2) +
  geom_pointrange() +
  geom_line() +
  scale_x_continuous(breaks = 2013:2024) +
  theme_minimal()

1,532 trillion

pop_mult = 1.008

df_temp_pop <- data.frame(
  year = c(2019, 2020, 2021, 2022, 2023, 2024, 2025, 2026, 2027, 2028, 2029, 2030, 2031, 2032, 2033, 2034, 2035),
  resident = c(4026209, 4044210, 3986842, 4073239, 4149253, 4180868, 4204515, 4204515 * pop_mult,
                 4204515 * pop_mult ** 2, 4204515 * pop_mult ** 3, 4204515 * pop_mult ** 4, 4204515 * pop_mult ** 5, 4204515 * pop_mult ** 6, 4204515 * pop_mult ** 7, 4204515 * pop_mult ** 8, 4204515 * pop_mult ** 9, 4204515 * pop_mult ** 10),
  population = c(5703569, 5685807, 5453566, 5637022, 5917648, 6036860, 6111175, 6111175 * pop_mult,
                 6111175 * pop_mult ** 2, 6111175 * pop_mult ** 3, 6111175 * pop_mult ** 4, 6111175 * pop_mult ** 5, 6111175 * pop_mult ** 6, 6111175 * pop_mult ** 7, 6111175 * pop_mult ** 8, 6111175 * pop_mult ** 9, 6111175 * pop_mult ** 10)
)

ggplot(df_temp_pop, aes(x = year, y = population)) +
  geom_line(alpha = 0.5, colour = "royalblue", size = 1.2) +
  geom_point(alpha = 0.5, colour = "royalblue", size = 3) +
  geom_line(data = df_temp_pop %>% filter(year <= 2025), colour = "royalblue", size = 1.2) +
  geom_point(data = df_temp_pop %>% filter(year <= 2025), colour = "royalblue", size = 3) +
  theme_minimal()

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.